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Abstract 

A variation of Gallager error-correcting codes is investigated using statistical 
mechanics. In codes of this type, a given message is encoded into a codeword 
which comprises Boolean sums of message bits selected by two randomly con- 
structed sparse matrices. The similarity of these codes to Ising spin systems 
with random interaction makes it possible to assess their typical performance 
by analytical methods developed in the study of disordered systems. The typ- 
ical case solutions obtained via the replica method are consistent with those 
obtained in simulations using belief propagation (BP) decoding. We discuss 
the practical implications of the results obtained and suggest a computation- 
ally efficient construction for one of the more practical configurations. 
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I. INTRODUCTION 



Error-correcting codes are commonly used for reliable data transmission through noisy 
media, especially in the case of memoryless communication where corrupted messages cannot 
be repeatedly sent. These techniques play an important role in a wide range of applications 
from memory devices to deep space explorations, and are expected to become even more 
important due to the rapid development in mobile phones and satellite-based communica- 
tion. 

In a general scenario, the sender encodes an dimensional Boolean message vector ^, 
where C,i € (0, 1), Vi, to an M(> A^) dimensional Boolean codeword zq, which is then being 
transmitted through a noisy communication channel. Noise corruption during transmission 
can be modelled by the noise vector ^, where corrupted bits are marked by the value 1 
and all other bits are zero, such that the received corrupted codeword takes the form z = 
Zq + C (mod 2). The received corrupted message is then decoded by the receiver for retrieving 
the original message ^. 

The error-correcting ability comes at the expense of information redundancy. Shannon 
showed in his seminal work |l| that error-free communication is theoretically possible if 
the code rate, representing the fraction of informative bits in the transmitted codeword, is 
below the channel capacity; in the case of unbiased messages transmitted through a Binary 
Symmetric Channel (BSC), which we will focus on here, R = N/M satisfies 



The expression on the right is termed Shannon's bound. However, Shannon's derivation is 
non-constructive and the quest for codes which saturate Eq.(|l|) has been one of the central 
topics of information theory ever since. 

In this paper we examine the efficiency and limitations of Gallager-type error-correcting 
code 0,0, which attracted much interest recently among researchers in this field. This code 
was discovered almost forty years ago by Gallager but was abandoned shortly after its 




(1) 
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invention due to the computational limitations of the time. Since their recent rediscovery 
by MacKay and Neal 0, different variations of Gallager-type codes have been developed 
attempting to get as close as possible to saturating Shannon's bound. 
In this paper we will examine the typical properties of a family of codes based on one 



variation, the MN code using the established methods of statistical physics P|-pi|], to 
provide a theoretical study based on the typical performance of codes rather on the worst 
case analysis. 

This paper is organised as follows: In the next two sections, we introduce Gallager-type 
error-correcting codes in detail and link them to the statistical mechanics framework. We 
then examine the equilibrium properties of various members of this family of codes using 
the replica method (section IV) and compare the bit-error rate below criticality. In section 
V, we examine the relation between Belief-Propagation (BP) decoding and the Thouless- 
Anderson-Palmer (TAP) approach to diluted spin systems; we then use it for comparing 
the analytical results obtained via the replica method to those obtained from simulations 
using BP decoding. In section VI we show a computationally efficient construction for one 
of the more practical constructions. Finally, we present conclusions for the current work and 
suggest future research directions. 



II. GALLAGER-TYPE ERROR-CORRECTING CODES 

There are several variations in Gallager-type error- correcting codes. The one discussed 
in this paper is termed the MA^code, recently introduced by MacKay and Neal 0. In these 
codes, a Boolean message ^ is encoded into a codeword Zq using two randomly constructed 
Boolean sparse matrices Cg and Cn, which are characterised in the following manner. 

The rectangular sparse matrix Cs is of dimensionality M x N, having randomly chosen 
K non-zero unit elements per row and C per column. The matrix C„ is an M x M (mod 2)- 
invertible matrix having randomly chosen L non-zero elements per row and column. These 
matrices are shared by the sender and the receiver. 
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Using these matrices, one can encode a message ^ into a codeword Zq in the following 
manner 

z, = C-'Cs^ (mod 2), (2) 

which is then transmitted via a noisy channel. Note that all matrix and vector components 
are Boolean (0,1), and all summations are carried out in this field. For simplicity, the 
noise process is modelled hereafter by a binary symmetric channel (BSC), where each bit is 
independently flipped with probability p. Extending the code presented here to other types 
of noise is straightforward. 

During transmission, a noise vector ^ is added to Zq and a corrupted codeword z = 2o + C 
(mod 2) is received at the other end of the channel. Decoding is then carried out by taking 
the product of the matrix Cn and the received codeword z, which results in Cs^ + CnC = 
CnZ = J. The equation 

CsS + CnT = J (mod 2), (3) 

is solved via the iterative methods of belief propagation (BP) [0,0 to obtain the most 
probable Boolean vectors S and r. BP methods in this context have recently been shown 



to be identical to a TAP |]I4| based solution of a similar physical system 



III. A STATISTICAL MECHANICS PERSPECTIVE 

Sourlas was the first to point out that error-correcting codes of this type have a similarity 
to Ising spin systems in statistical physics he demonstrated this using a simple version 



of the same nature. His work, that focused on extensively connected systems, was recently 
extended to finitely connected systems [^|11|. We follow a similar approach in the current 
investigation; preliminary results have already been presented in |10| . 



To facilitate the statistical physics analysis, we first employ the binary representation 
(±1) of the dynamical variables S and r and of the check vector J rather than the Boolean 
one (0, 1). The /i-th component of Eq.(^ is then rewritten as 



n ^ n 



J. 



(4) 



where Cs{fi) and Cn{fJ') are the sets of all indices of non-zero elements in row fi of the 
sparse matrices Cg and C„, respectively. The check /i is given by message ^ and noise C as 
jeCnM Oi should be emphasised that the message vector ^ and the noise 
vector ^ themselves are not known to the receiver. 

An interesting link can now be formulated between the Bayesian framework of MN 
codes and Ising spin systems. Rewriting Kronecker's delta for binary variables x and y 
as = (1 + xy)/2 = lim^_»oo exp(— 1; a;?/]), one may argue that, using it as a 

likelihood, equation gives rise to the conditional probability of the check J for given S, 
r, C, and Cn 

M 



P(J|5,T,C„C„)= limexp -/3 5]5[-l;J,, 11 11 ^j] 



(5) 



Prior knowledge about possibly biased message and noise is represented by the prior distri- 
butions 



VJS) 



exp {FsEIiS, 



(2 cosh F,)^ ' ' (2 cosh F, 

respectively. Non-zero field Fg is introduced for biased message and F„ is determined by 
flip rate p of channel noise as F„ = (1/2) In ((1 ~p)/p)- Using equations (|) and (^), the 
posterior distribution of S and r for given check J and matrices Cs and C„ is of the form 

ns,r\j,c^,c.) = riJ\c.,c.) 

^ ^.^ exp{-(3n{S,T\J,V)) 

Z{J,V) 

where P( J|C„ C^) = E^S,t} ^(-^I"^' ^' a)n(S)P„(T), 

Jk 

— 1; >J<h,..,iK;ji,--,jL 



Vn{T) 



exp 



(6) 



(7) 



7^(5, r\J,V) = Y.5[-l■,J,X{S,\{r,]-^Y.S^-^Y: r, 



~ $Z '^<h,--,iK;ji,-,jL>^ 
<h,--,iK;ji,--,jL> 

p N p M 



>Sii ■ ■ ■ ^iK'^ji ■ ■ ■ '^j 



and 



Z{J,V)= lim exp{-/3n{S,T\J,V)) 



^ "I" rJ^{h,---,iK\ii,---,3L) i'^<h,--,iK;ji,--,jL>^h • • • ^iii'^ji • • • 'Ol -'-) 

-> 



= E n 

{S,t} <«i,--,«K;ii,-- Jl'-; 

/ N M 

xexp F,5]^, + F„^r, I . (9) 



The final form of posterior distribution (j^) implies that the MA'^code is identical to an Ising 
spin system defined by the Hamiltonian (P) in the zero temperature limit T = f3~^ 0. 
In equations (Q) and (j^), we introduced the sparse connectivity tensor ^^<ji,..,j^> which 
takes the value 1 if the corresponding indices of both message and noise are chosen (i.e., if 
all corresponding indices of the matrices Cs and Cn are 1) and otherwise, and coupling 
,iK;ji,--,jL> ~ ^h^h ■ ■ -^iifOiCj^ • • • Cjl- These come to isolate the disorder in choosing 
the matrix connections, embedded in 1^<ii,..,jL>, simplify the notation. 

The posterior distribution (0) can be used for decoding. One can show that expectation 
of the overlap between original message ^ and retrieved one ^ 

1 ^ ^ 

^ = i7E^^e., (10) 

i=l 



is maximised by setting ^ to its Bayes-optimal estimator |T6|-19| 



ef = sign(mf), mf= Y: 5. r| J, C., a). (11) 

{S,T} 

It is worth while noting that this optimal decoding is realized at zero temperature rather 
than at a finite temperature as in |[T^-|19|. The reason is that the true likelihood term 
corresponds to the ground state of the first term of the Hamiltonian due to the existence 
of more degrees of freedom, in the form of the dynamical variables r, which do not appear 
in other systems. Introducing the additional variables r, the degrees of freedom in the 
spin system increase from to + M, while the number of constraints from the checks 
J remains M. This implies that in spite of the existence of quenched disorder caused by 
JT" and V, the system is free from frustration even in the low temperature limit, which is 
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useful for practical decoding using local search algorithms. The last two terms in Eq.(^ 
scale with (3 remain finite even in the zero temperature limit (3 —>■ oo representing the true 
prior distributions, which dominates the statistical properties of the system, while the first 
term vanishes to satisfy the parity check condition (^. 

IV. EQUILIBRIUM PROPERTIES: THE REPLICA METHOD 

As we use the methods of statistical mechanics, we concentrate on the case of long 
messages, in the of limit of A^, M — oo while keeping code rate R = N/M = K/C finite. 
This limit is quite reasonable for this particular problem since Gallager-type codes are usually 
used in the transmission of long (10^ — 10^) messages, where finite size corrections are likely 
to be negligible. 

Since the first part of the Hamiltonian (^) is invariant under the gauge transformation 
Si — ^ iiSi, Tj (jTj and 1, it is useful to decouple the correlation between 

the vectors S, r and ^, C- Rewriting the Hamiltonian using this gauge, one obtains a 
similar expression to Eq.(|I) apart from the second terms which become Fg/ f3J2i=i^i'Si and 

Fn/PJ2j=lCjTj- 

Due to the existence of several types of quenched disorder in the system, it is natural 
to resort to replica method for investigating the typical properties in equilibrium. More 
specifically, we calculate expectation values of n-th power of partition function with 
respect to the quenched variables ^, ^ and V and take the limit n — 0. 

Carrying out the calculation in the zero temperature limit (3 oo gives rise to a set of 
order parameters 

qa,p,..,, = (^j:Z.SrS^,..,Sj) , r„,,,..,,= (-^f;y,r;r/,..,r7\ (12) 

where a, (3, .. represent replica indices, and the variables Zi and Yj come from enforcing the 
restriction of C and L connections per index, respectively ■ 

YP<..'. ..> - c] = f zS<.......f , (13) 

,. . \ Jo Zn 



and similarly for the restriction on the j indices. 

To proceed further, it is necessary to make an assumption about the order parame- 
ters symmetry. The assumption made here is that of replica symmetry in both the order 
parameters and the related conjugate variables 

9a,/3..7 ~ ^1 J ^(^) ' ^"./9--7 ~ '^g J ^(S:) (14) 
ra,p..-y = ar J dy p{y) , ra,p..^ = a- J dy p{y) y^ , 

where / is the number of replica indices, are normalisation coefficients, and vr(x), tt (x), p{y) 
and p{y) represent probability distributions. Unspecified integrals are over the range 
[—1,+!]. This ansatz is supported by the facts that (i) the current system is free of frus- 
tration and (ii) there has never been observed replica symmetry breaking at Nishimori's 
condition 



which corresponds to using correct priors and F„ in our case [|T6|. The 
results obtained hereafter also support this ansatz. Extremizing the partition function with 
respect to distributions 7r(-), 7r(-), p(-) and p(-), one then obtains the free energy per spin 

( C r CL r 

= extr{^,^-,p,p} l—ln2 + C J dx dx it{x) 7r(x) ln(l + xx) + — J dy dy p{y) p{y) ln(l + yy) 



K 



K 



W dXk1T{Xk) 
k=l 

n dxkTf{xk 



c 

n 

.k=l 



c_ 

K 



c 



n dyipim 



.1=1 



L 1 r K L - 

\[dyip{,yi) In i+n^^n^/' 

1=1 J L k=l 1=1 . 

In e^^^f{{l + Xk)+e-^^^\{{l~x,) 

k=l k=l 

In e^"^n(l + y,) + e-^"^n(l-^/) 
/=i 1=1 



(15) 



where angled brackets with subscript ^, <^ and T) denote averages over the message and noise 
distributions respectively, and sparse connectivity tensor V. Message averages take the form 

1 + ^ tanh Fs 



(■■■),= E 



(16) 



€=±1 



and similarly for (■ ■ ■)^. Details of the derivation are given in Appendix A. 

Taking the functional variation of / with respect to the distributions vr, vr, p and p, one 
obtains the following saddle point equations 
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, C-1 



C-1 



K-l L / K-l L \ 

n dxi 7r(xi) / n dyi p{yi)5 [x - H 11 ' 



p{y) 
piy) 



1=1 

L-l 



n dyi piyi)5 
1=1 \ 

n (^yi Piyi) (^^ (^y ~ ^^^^ (^CFn + J2 tanh^^ yi 

.L-l I K L-l \ 

J n (^yi p(yi)^ [y ~I[^iI[yn ■ 



1=1 

K 

Wdxi tt{xi 
1=1 



1=1 1=1 

L-l 



(17) 



1=1 



1=1 1=1 



After solving these equations, the expectation of the overlap between the message $, and the 
Bayesian optimal estimator (pT]), which serves as a performance measure, can be evaluated 

as 

1 / \ 



T?(E^*sign {S,)^^) 



dz (f){z) sign(z) 



(18) 



where 



c 



Y[dxi Tf{xi) 
ll=l 



6 z 



tanh (f,^ + tanh"^ 
V 1=1 



(19) 



The derivation of Eqs.(rB) and (|T^) is given in Appendix B. 

Examining the physical properties of the solutions for various connectivity values exposes 
significant differences between the various cases. In particular, these solutions fall into three 
different categories: the cases of i^' = 1 and general L value, the case of = L = 2 and all 
other parameter values where either > 3 or L > 3 (and > 1). We describe the results 
obtained for each one of these cases separately. 



A. Analytical solution - the case of K > 3 or L > 3, K > 1 

Results for the cases of-ft'>3orL>3,-ft'>l can be obtained analytically and have 
a simple and transparent interpretation; we will therefore focus first on this simple case. 
For unbiased messages (with Fg = 0), one can easily verify that the ferromagnetic phase, 
characterised by m = 1, and the probability distributions 

7r(x) = 6{x - 1), 7r(a;) = 6{x - 1), p{y) = 6{y - 1), p{y) = 6{y - 1) ; (20) 
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and the paramagnetic state of m = with the probabihty distributions 

tt{x) = 5{x), Tf{x) = 5{x), piy) = 5{y), 

, , l + tanhF„^, , , 1— tanhF„^, , , , , 

P{V) = ^ -^iy - tanhF„) + ^6{y + tanhF„), (21) 



satisfy saddle point equations (p!7|) . Other solutions may be obtained numerically; here we 
have represented the distributions by 10'^ — 10^ bins and iterated Eqs.([T7|) 100 — 500 times 
with 10^ Monte Carlo sampling steps for each iteration. No solutions other than the above 
two have been discovered. 

The thermodynamically-dominant state is found by evaluating the free energy of the two 
solutions using Eq.(|T^), which yields 

C 1 
/ferro = - "^^n tanh F„ = - — F„ tanh F„ , (22) 

for the ferromagnetic solution and 

/para = ^ In 2 - In 2 - ^ In 2 cosh F„ = In 2 - In 2 - ^ In 2 cosh (23) 

for the paramagnetic solution. 

Figure 0(a) describes schematically the nature of the solutions for this case, in terms 
of the free energy and the magnetisation obtained, for various flip rate probabilities. The 
difference between the free energies of Eqs. (^2D and ( p3|) 



/ferro " /para =^^[R-l + H,{p)] , (24) 

vanishes in the boundary between the two phase 

R, = l-H2{p) = l +p\og,{p) + {l-p) \og,{l - p), (25) 

which coincides with Shannon's channel capacity. 

Equation ( P^D indicates that all constructions with either > 3 or L > 3 (and K > 1) 
can potentially realize error-free data transmission for R < Rc in the limit where both 
message and codeword lengths and M become infinite, thus saturating Shannon's bound. 
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B. The case of K = L = 2 



All codes with either K = 3otL = 3, K>1 potentially saturate Shannon's bound and 
are characterised by a first order phase transition between the ferromagnetic and paramag- 
netic solutions. On the other hand, numerical investigation based on Monte Carlo methods 
indicates of significantly different physical characteristics for K = L = 2 codes shown in 
Fig.|I|(b). 

At the highest noise level, the paramagnetic solution (^Tj) gives the unique extremum of 
the free energy until noise level reaches the first critical point pi, at which the ferromagnetic 
solution (^) of higher free energy appears to be locally stable. As the noise level decreases, 
a second critical point p2 appears, where the paramagnetic solution becomes unstable and a 
sub-optimal ferromagnetic solution and its mirror image emerge. Those solutions have lower 
free energy than the ferromagnetic solution until the noise level reaches the third critical 
point J93. Below p^, the ferromagnetic solution becomes the global minimum of the free 
energy, while the sub-optimal ferromagnetic solutions still remain locally stable. However, 
the sub-optimal solutions disappear at the spinodal point ps and the ferromagnetic solution 
(and its mirror image) becomes the unique stable solution of the saddle point Eqs.(|l^ as 
shown by the numerical investigation for all p < Ps- 

The analysis implies that p^, the critical noise level below which the ferromagnetic solu- 
tion becomes thermodynamically dominant, is lower thanpc = H2^{l—R) which corresponds 
to Shannon's bound. Namely, K = L = 2 does not saturate Shannon's bound in contrast 
to K > 3 codes even if optimally decoded. Nevertheless, it turns out that the free energy 
landscape, for noise levels < p < Ps, offers significant advantages in the decoding dynamic 
comparing to that of other codes {K >3otL>3,K>1). 
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C. General L codes with K = 1 



The particular choice of = 1, independently of the value chosen for L, exhibits a 
different behaviour presented schematically in Fig.|l](c); also in this case there are no simple 
analytical solutions and all solutions in this scenario, except for the ferromagnetic solution, 
have been obtained numerically. The first important difference to be noted is that the 
paramagnetic state (^Tj) is no longer a solution of the saddle point equations (P!7| ) and is 
being replaced by a sub-optimal ferromagnetic state. Convergence to the perfect solution of 
m = 1 can only be guaranteed for corruption rates smaller than that of the spinodal point, 
marking the maximal noise level for which only the ferromagnetic solution exists, p < ps. 

The K = 1 codes do not saturate Shannon's bound in general; however, we have found 
that at rates i? < 1/3 they outperform the K = L = 2 code while offering slightly improved 
dynamical (decoding) properties. Studying the free energy in this case shows that as the 
corruption rate increases, sub-optimal ferromagnetic solutions (stable and unstable) emerge 
at the spinodal point ps- When the noise increases further this sub-optimal state becomes 
the global minimum at pi, dominating the system's thermodynamics. The transition at pi 
must occur at noise levels lower or equal to the value predicted by Shannon's bound. In 
Fig.0 we show free energy values computed for a given code rate and several values of L, 
marking Shannon's bound by a dashed line; it is clear that the thermodynamical transition 
observed numerically (i.e. the point where the ferromagnetic free energy equals the sub- 
optimal ferromagnetic free energy) is bellow, but very close, to the channel capacity. It 
implies that these codes also do not quite saturate Shannon's bound if optimally decoded 
but get quite close to it. 

V. DECODING: BELIEF PROPAGATION/TAP APPROACH 

The Bayesian message estimate (^^) potentially provides the optimal retrieval of the 
original messages. However, it is computationally difficult to follow the prescription exactly 
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as it requires a sum over 0{2^) terms. Belief propagation [ll2| , p!3| (BP) can be used for 
obtaining an approximate estimate. It was recently shown |p that the BP algorithm can be 
derived, at least in the current context, from the TAP approach [|T4| to diluted systems in 
statistical mechanics. 

Both algorithms (BP/TAP) are iterative methods which effectively calculate the 
marginal posterior probabilities V{Si\J,Cs,Cn) = I]{{5fc^j,T} "^(-^j ''"l-^) C's, C„) and 

{Tfe^j}} '^i^^ 'T\J 1 Cs, Cn) based on the following three assumptions: 

1. The posterior distribution is factorizable with respect to dynamical variables Si=i^,,,^N 
and rj=i^...^M- 

2. The influence of check J^=i,...,Af on a specific site Si (or Tj) is also factorizable. 

3. The contribution of a single variables 5*4=1^. ..^at, rj=i^..,^M and J^=i^,,,^m to the macro- 
scopic variables is small and can be isolated. 

Parameterising pseudo-marginal posteriors and marginalized conditional probabilities as 

^(5,,|{J.^^},C„C„) = i^|^^, p(r,|{J.^^},C„C„) = i^|^, (26) 

V{J,\S,, { J.^^}, C., Cn) ~ i±^^^, P( J^|r„ { J.^^}, C., Cn) ~ i±|^, (27) 

the above assumptions provide a set of self-consistent equations P,pi| 

mil = tanh ( + ^ tanh~^(m^;) ) , m^J; = tanh ( F„ + ^ tanh-^(m;;;) ) . (2^ 



and 



</ = '^M n n ^'"^i = jf^ n n 



Here, Ais{l) and A^„(/) indicate the set of all indices of non-zero components in the Z-th 
column of the sparse matrices Cg and C„, respectively. Similarly, Cs{fi) and Cn{lj) denote 
the set of all indices of non-zero components in /x-th row of the sparse matrices Cg and Cn, 
respectively. The notation Cs{fi)/l represents the set of all indices belonging to i2s(/i) except 
the index /. 
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Equations (pH]) and (pQ]) are solved iteratively using the appropriate initial conditions. 
After obtaining a solution to all m^; and m^;, an approximated posterior mean can be 
calculated as 



tanh F5 + ^ tanh ""^i 



m 



fit J 



(30) 



which provides an approximation to the Bayes-optimal estimator ( pH]) in the form of = 
sign(mf ). 

Notice that the rather vague meaning of the fields distributions introduced in the previous 



section becomes clear by introducing the new variables x = C,im 



fa ' 



X 



and y = Cj^^j iHl- If one considers that these variables are independently drawn from the 
distributions vr(a;), n{x), p{y) and p(y), the replica symmetric saddle point equations (|T^ 
are recovered from the BP/TAP equations ( p8|) and (^9|) . This connection can be extended 
to the free energy as equations (^) and (^91) extremize the TAP free energy 

^ 1 A/ , M 



A/f 1 1 

/TAp({m},{rh}) = -ln2 + -E E 1- (l + + ]^ E E 1- (l + 



M 



-i^Ein n < n 

1 ^ 
--Vln 



m, 



M 



-Vln 

AT ^ 



N 



n 

fieMsii) 

n 



+e- 



n 

Ai6A1s(i) 



n (1 

AtGA1„(j) 



(31) 



This expression may be used for selecting the thermodynamically dominant state when 
Eqs.([28|) and (pOj) have several solutions. 

We have investigated the performance of the various codes using BP/TAP equations as 



the decoding algorithm. Solutions have been obtained by iterating the equations ( pSj ) and 
(p9|) 100 — 500 times under various initial conditions. Since the system is not frustrated, the 
dynamics converges within 10 — 30 updates in most cases except close to criticality. The 
numerical results mirror the behaviour predicted by the analytical solutions. 
For either i^>3orL>3,-ft'>l codes, the ferromagnetic solution 
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= 6, Ki = Kj = 0, Ki = 0, (32) 

which provides perfect decoding (m = 1) and the paramagnetic solution (m = 0) 

m% = Q, m% = Q, m;:^. =tanhF„ = l-2p, = Q , (33) 

are obtained in various runs depending on the initial conditions (the message is assumed 
unbiased resulting in = 0). However, it is difficult to set the initial conditions within 
the basin of attraction of the ferromagnetic solution without prior knowledge about the 
transmitted message 

Biased coding is sometimes used for alleviating this difficulty 0. Using a redundant 
source of information (equivalent to the introduction of a non-zero field Fs in the statis- 
tical physics description), one effectively increases the probability of the initial conditions 
being closer to the ferromagnetic solution. The main drawback of this method is that the 
information per transmitted bit is significantly reduced due to this redundancy. In order to 
investigate how the maximum performance is affected by transmitting biased messages, we 
have evaluated the critical information rate (i.e., code rate xH2{fs = (1 + tanhF5)/2), the 
source redundancy), below which the ferromagnetic solution becomes thermodynamically 
dominant [Fig.^(a)]. The data were obtained by the BP/TAP method (diamonds) and nu- 
merical solutions of from replica framework (square); the dominant solution in the BP/TAP 
results, was selected by using the free energy (0). Numerical solutions have been obtained 
using 10^ — 10"^ bin models for each distribution and had been run for 10^ steps per noise 
level. The various results are highly consistent and practically saturate Shannon's bound 
for the same noise level. However, it is important to point out that close to Shannon's 
limit, prior knowledge on the original message is required for setting up appropriate initial 
conditions that ensure convergence to the ferromagnetic solution; such prior knowledge is 
not available in practice. 

Although K, L > 3 codes seem to offer optimal performance when highly biased messages 
are transmitted, this seems to be of little relevance in most cases, characterised by the 
transmission of compressed unbiased messages or only slightly biased messages. In this 
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sense, K = L = 2 and K = 1 codes can be considered more practical as the BP/TAP 
dynamics of these codes exhibit unique convergence to the ferromagnetic solution (or mirror 
image in the K = L = 2 case) from any initial condition up to a certain noise level. This 
property results from the fact that the corresponding free energies have no local minima 
other than the ferromagnetic solution below Ps- 

In figures ^(b) and (c) we show the value of ps for the cases oi K = L = 2 and K = 1, 
L = 2 respectively, evaluated by numerical solutions from the replica framework (diamonds) 
and using the BP/TAP method. 

The case of K = L = 2 shows consistent successful decoding for the code rates examined 
and up to noise levels slightly below, but close to. Shannon's bound. It should be emphasised 
here that initial conditions are chosen almost randomly in the BP/TAP method, with a very 
slight bias of (9(10^^^) in the initial magnetisation. This result suggests using K = L = 2 
codes (or similar), rather than K,L > 3 codes, although the latter may potentially have 
better equilibrium properties. 

In Fig.|^(c) we show that for code rates R < 1/3, codes parametrised by = 1 and 
L = 2 outperform K = L = 2 codes with one additional advantage: Due to the absence of 
mirror symmetries these codes converge to the ferromagnetic state much faster, and there is 
no risk of convergence to the mirror solution. The difference in performance becomes even 
larger as the code rate decreases. Higher code rates will result in performance deterioration 
due to the low connectivity, eventually bringing the system below the percolation threshold. 

In Fig.^ we examine the dependence of the noise level of the spinodal point Ps on the 
value of L, and show that the choice of L = 2 is optimal within this family. Codes with 
L = 1 have very poor error-correction capabilities as their Hamiltonian (|]) corresponds to 
the Mattis model, which is equivalent to a simple ferromagnet in a random field attaining 
magnetisation m = 1 only in the noiseless case. 
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VI. REDUCING ENCODING COSTS 



The BP/TAP algorithm already offers an efficient decoding method, which requires 0{N) 
operations; however, the current encoding scheme includes three costly processes: (a) The 
computational cost of constructing the generating matrix C^^Cs requires 0{N^) operations 
for inverting the matrix C„ and 0{N'^) operations for the matrix multiplication, (b) The 
memory allocation for generating the matrix C^^Cs scales as 0{N'^) since this matrix is 
typically dense, (c) The encoding itself Zq = C^^Cs$, (mod 2) requires 0(A^^) operations. 

These computational costs become significant when long messages = 10^ ~ 10^ are 
transmitted, which is typically the case for which Gallager-type codes are being used. This 
may require long encoding times and may delay the transmission. 

These problems may be solved by utilising systematically constructed matrices instead 
of random ones, of some similarity to the constructions of Here, we present a simple 
method to reduce the computational and memory costs to 0{N) for K = L = 2 and K = 1, 
L = 2 codes. Our proposal is mainly based on using a specific matrix for C„, 
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instead of a randomly-constructed one. For Cg, we use a random matrix of K = 2 {ot K = 1) 
non-zero elements per row as before. 

The inverse (mod 2) of becomes the lower triangular matrix 
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10 
110 
1110 

1111 



■■0 
■■0 
■■0 
■■0 



1111 



1 1 



(35) 



This suggests that encoding the message $, into a codeword z° would require only 0{N) 
operations by carrying it out in two steps 



t^ = (C.OM (mod 2), for/i = l,2,...,M, 

4 = iC-H)^ = z%+t, (mod 2), for/i = 2,...,M, 



(36) 
(37) 



with z1 = ti. Both steps require 0{N) operations due to the sparse nature of C^. In 
addition, the required memory resources are also reduced to 0{N) since only the sparse 
matrix Cs should be stored. 

The possible drawback of using the systematic matrix (|34D is a deterioration in the error 
correction ability. We have examined numerically the performance of new construction to 
discover, to our surprise, that it is very similar to that of random matrix based codes as 
shown in Table |. Although our examination is only limited to BSC and i.i.d. messages, 
it seems to suggest that some deterministically constructed matrices may be implemented 
successfully in practice. 



VII. SUMMARY 

In this paper, we have investigated the typical performance of the MN codes, a variation 
of Gallager-type error- correcting codes, by mapping them onto Ising spin models and making 
use of the established methods of statistical physics. We have discovered that for a certain 
choice of parameters, either A">3orL>3,A'> 1 these codes potentially saturate the 
channel capacity, although this cannot be used efficiently in practice due to the decrease in 
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the basin of attraction which typically diverts the decoding dynamics towards the undesired 
paramagnetic solution. 

Codes with K = 2 and L = 2 show close to optimal performance while keeping a large 
basin of attraction, resulting in more practical codes. Constructions of the form K — 1, 
L — 2 outperform the K — L — 2 codes for code rates it! < 1/3, having improved dynamical 
properties. 

These results are complementary to those obtained so far by the information theory 
community and seem to indicate that worst-case analysis can be, in some situations, too 
pessimistic when compared to the typical performance results. 

Beyond the theoretical aspects, we proposed an efficient method for reducing the com- 
putational costs and the required memory allocation by using a specific construction of the 
matrix C„. These codes are highly attractive and provide lower computational costs for 
both encoding and decoding. 

Various aspects that remain to be studied include a proper analysis of the finite size effects 
for rates below and above the channel capacity, which are of great practical relevance; and 
the use of statistical physics methods for optimising the matrix constructions. 
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APPENDIX A: REPLICA FREE ENERGY 



The purpose of this appendix is to derive the averaged free energy per spin (p^f). Applying 
the gauge transformation 



(Al) 



to eq. @, one may rewrite the partition function in the form 



N M 

Zit CV) = Y: exp I Fs ^ ^.Si + Fr,J2C: 



i=l 



X n 



(A2) 



Using the rephca method, one calculates the quenched average of the n-th power of the 
partition function given by 



xi-- x" \ V i=l a=l / I ^\ \ j=l a=l 

n n 1 + ^^^(n,-....-..) (^n ■ ■ ■ ^1^^ ■ ■ ■ ^ - 



V 



where averages with respect to £ can be easily performed 

/ (^^^^\\ JI r/l + tanhF,\ p v" /l-tanhF,\ p v" 
(exp(F,E6E^r)) =n[(^ 2 ^je^^S...^. ije-^^S. 



(A3) 



i=l a=l 



COL 

a-1 "^i 



= n (^exp (^ei^. E j ^ , (A4) 

and similarly for (■ ■ ■)^. The main problem is in averages over the sparse tensor realisations 
which have complicated constraints. Following the procedure introduced by Wong and 
Sherrington it is being rewritten as 



n n 1 + 2^{iu-,iK;h,-,jL) {si ■ ■ ■ s^^Tf^ ■ ■ ■ - 1 

^(*i,---,ix;ii,---,i_L> "=i 
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■^'^ HW^y ^(i,i2,-,ix;ji,-,ii,> - n ^ I ^(ii,-,iK;ij2,-jL> 

V 1=1 \{ii,---,iK;ji,---,jL) I i=i ;ii,---,jL) 

n r 



X n n 



L 



(A5) 



where ■ ■) represents Dirac's 5-function and 

= E n 4 E ^^(..2,-.x;.. -..> - ^ ) n 4 ^ ^^(n,-.x;..2, -,,.> - ^ | (A6) 

V i-\ \{i2,-,iK]ji,-,jL) I i=l \{ii,-,iK]j2,-,jL) 

represents the normahsation constant. 

We first evaluate this normalisation constant using the integral representation of the 
5-function and Eg. ([A6|) , to obtain 



N 



V i=l \{i2,-,iK]jl,-,jL) 



M I 

I n ^ I E ^(n,-,iif;i,i2, -,ii,> " L 



(«2,---,«K;ii,---,jL> 



n / ^e--Mn / B 



^ [jo 2ti 



3 ^—iLvn 



N 

xEn 



M 



D j=l \^{i2,-,iK]jl,---,jL) ) J=l L (*l''--'«K;j2,---,iL) 



^4 r z-^- dz/, 



n / ^e--Mn / B 



l^io 27r 



J —iLuj 



1=1 V— J J 

xE n {(e"' 
^ («i,---,«A';ii,---jL> 



1 . . . e*^«K e*''^! • • • e'^'L 



27ri ^' 



n (l + Z,,---Z,,F,,---r,J, (A7) 
(«i,---,«A';ii,---,iL> 



where we made use of the transformations Zi = e ' , Yj = e^^^ , and carried out summations 
with respect to the realisation of T). Expanding the product on the right hand side one 
obtains 



(nr--:«A;ii>--- jl> 
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exp 



exp 



exp 



Y: \n{l + {Z,,...Z,,Y,,...Y,,)} 

_{h,---,iK;ji,---,jL) 



.(«i.---.«x;iir--jL> 



3Lt 



N 



K 



\i=l 



1 

L! 



M 



in the thermodynamic hmit. Using the identities 



1 = Jdqs(j2Z^-q], 1= J dr 6(^2^3 

\i=l ) \j=l 



Eq. ( |A7D becomes 

M= Jdq6(j2Zi-(l) Jdrslj^yj 

\i=l / \7=1 



N 

xn 

1=1 



dZi „_ 



£^-(C+l) 



27ri 
dq 



M 

n 



27ri ^ 



exp 



2m 



exp 



TV 



(if 



AT 



27ri * 



M 

n 



27ri 
2-Ki ^ I 



exp 



exp 



M 



27ri 



dr 



K T 



dq I : I dr I : exp — — — qq — rf 



27ri 



N 

xn 

1=1 



—Z^' exp [qZ, 



K\ L\ 

M 

n 

J i=l L' 



exp {fYA 

2'Ki ■' 



The contour integrals provide the following constants 



(A8) 



(A9) 



N 

n 

i=l 



(C+1) ^^^^ 

Ztxi 



exp (fF,) 



27ri 



(AlO) 



respectively. Applying the saddle point method to the remaining integrals, one obtains 



M = extT{q^q^r,f} { GXp 



qq-rf + NClnq-N ln(C!) + MLlnf - M ln{L\] 



(A12) 



which yields the following saddle point equations with respect to q, r, q and f 

NC ML 
q = — , r = — 



{k-i)\1j: 



{L-iy.m 



(A13) 
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(AM) 



providing the normalisation constant 

Equation ( |A5|) can be evaluated similarly. Following a similar calculation to that of 
Eq.(^) provides 



U=\^\ (^-J exp(^— --gg-rr^ 



3l 



X n 

{h,--,iK;ji,-,jL) 



a=l 



Using the expansion 

n 



■ T 



a=l 



a a 
iK ji 



a=l 



(ai,02> 



02 ^ 
ii > 



^"1 



V 



(A15) 



T- ■ ■ ■ T 



E E (^r 



resulting in 

n 



str) ■ ■ ■ {St:, ■ ■ ■ st-){T: 



°i . . . -7-"" 



f-r"! . . . T*^™ 



n 1 

l + (z,,•■■z,,r,,■■■r,J^^(l + ^^■■^^,r^ 



" . . . 



(A16) 



3 2" S 



(A17) 



Using the identities 



1 = J dqa,,..,a,^6 (^S', 

1 = J dro,,,..,aJ ^E^j 



St'-Z.-q 



(A18) 
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and going through the same steps as in Eqs. (|A9| - [A12|) , we arrive at 

r "1 

{il,---,iK;ji,---jL) L a-l 

= n n / dq.u-,e.J (E Sr--- St-Z, - q^r,..A 

m=0 {ai,---,am) \i=l I 

X / dT^,,..^^J ( f: • ■ ■ r;™r, - r«i,.,„ J exp ( ^ ( £ E " 



q.a \,- -,am ^Qi, -,( 



extr{q,q,r,f.}< exp 



On I E/ E/ 



' m=0 (ai,---,am> 



E/ E/ ''"ai,---,Om'^ai, 



-E E 

,---,a 

n N n 

+ E E ^oi,---,Qm E ^ ■ ■ ■ + E E ^ ■ J 

m=0 {ai,---,a„i) i=l m=0 (oi, ■■■,«„,) j=l 



M 



>. (A19) 



In order to proceed further, one has to make an assumption about the order parameter 
symmetry. We adopt here the rephca symmetric ansatz for the order parameters q, r, q and 
f. This imphes that the order parameters do not depend on the exphcit indices but only on 
their number. It is therefore convenient to represent them as moments of random variables 
defined over the interval [—1,1] 



qau-,»i = q J dx tt{x) x\ ra-^,-,ai = r J dy p{y) y\ 
gai,-,Q, = q dx 7c{x) x\ fa-^,-,ai = f dj) p{y) y\ 



(A20) 



Then, each term in Eg. ( |A19| ) takes the form 



E E '^^r^ '"^v r'" = E f!) / n dx,^{^^)x^ [ n dyMvr 



=0 m=0 k=l 



1=1 



E E 



K . K L / K L \" 

^— y n ^^fc 7r(x,,) j n dyi p{yi) 1 + H 11 1/' (^21) 

k=l 1=1 \ k=l 1=1 / 

gg E ( \ I dx dx 7r(x) ti{x)x^ x^ 

™ n \TfL / J 



qq / dx dx 7r(x) 7r(x) (1 + xx)" 



(A22) 



n N N . n 

E E qc.,,-,o.^Y.sr---srz, = qY.z,jdx7^{x)Y.s^'^ E ^r---^; 

Ti=0 i=l i=l m=0 (ai,---,Om) 



=0 (ai,---,Om) 
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N . n 

qY^Z, j dx7c{x) n (l + '^i 

i=l a=l 



X) 



Substituting these into ([A19|), one obtains 



(A23) 
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a=l ^ I C, 
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xextrj^^i^^p^^J exp 



27ri ^ / 



K L ^ 



— gg / dx dx ^{x) 7r(x)(l + xxY — rf dy dy p{y) p{y){l + yyY 



+qJ2Z^J dx T^ix) n (1 + Srx) +fY^Yjdy p{y) J] (1 + ^/y) 

1=1 a=l j=l a=l 

The term involving the spin variables 5* is easily evaluated using the residue theorem 



(A24) 



X exp 



l...Sni=i \ \ a=l //^j=lL"' 

fiC . C I n f C C 

yr, nldii ( n + £1) + 11(1 - X,) 

•' 1=1 \q=1 I 1=1 1=1 



i=l 



a=l 



/II, 
dx 7r{x) l[{l + S^x 

(A25) 



and similarly for the term involving the variables r. Substituting these into Eq. ( |A24|) , one 
obtains the n-th moment of partition function 



extrj^^^^p^^J exp 



-NC I^J dxdx 7r{x) 7r{x) ln(l + xx)"" - l| 



-ML ydydy p{y) p{y) ln(l + yyY - 1 
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n dyipivi) 
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(A26) 
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Finally, in the limit n — one obtains 
-(ln^(tC,P))«..= lto ,!n 



+ 



c_ 

K 



n (^yipivi 
.1=1 



In 



extr{^,ii-,p,p} |-— ln2 - C J dx dx n{x)7r (x) ln(l + xx) - — J dy dy p{y) p{y) ln(l + yij) 
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APPENDIX B: EVALUATION OF THE MAGNETISATION 



Here, we derive explicitly Eqs.(|T8D and (|T9|). After using the gauge transformation Si 
C,iSi, the magnetisation can be written as 



N 



m 



N 



^(sign(mi))|^^ 



fBll 



1=1 



introducing the notation nii = (Sj)^^^ (gauged average). 

For an arbitrary natural number p, one can compute p-th moment of rrii 



V 



lim lim 



E 



(B2) 



^{5\ri},...,{5",T"} ' 
where Tia denotes the gauged Hamiltonian of the a-th replica. Decoupling the dynamical 
variables and introducing auxiliary functions 7r(-), 7r(-), p(-) and /3(-), of a similar form to 
Eq. (|A20|), one obtains 



(B3) 



using the saddle point solution of 7r(-). 
Employing the identity 



n 



signfa;) = — 1 + 2 lim / 

n— >oo \ rn 

m=0 ^ 



2n\ /I + X 



2n—m 



1 — X 



(B4) 
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which holds for any arbitrary real number x G [—1, 1] and Eqs.( |B3| ) and ( |B4|) one obtains 



(sign(mi))|^^^^ =-1 + 2 j dz (j){z) jim ^ 
dz (f){z) sign(2;), 



m=0 



2n\ n+z 



m 



2n—m 



1 - Z 



(B5) 



where we introduced a new notation for the distribution 

c I c 



k=l 



4>{z) = Yldxi 7i:{xi) (6{z - - tanh ^ Xk) 
•' 1=1 \ 

thus reproducing Eqs.(|TB|) and (^). 



(B6) 
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FIGURES 




FIG. 1. Left hand figures show a schematic representation of the free energy landscape while 
figures on the right show the ferromagnetic, sub-optimal ferromagnetic and paramagnetic solutions 
as functions of the noise rate p; thick and thin lines denote stable solutions of lower and higher free 
energies respectively, dashed lines correspond to unstable solutions. {&) K > 3 or L > S, K > 1; the 
solid line in the horizontal axis represents the phase where the ferromagnetic solution (F, m = 1) 
is thermodynamically dominant, while the paramagnetic solution (P, m = 0) becomes dominant 
for the other phase (dashed line). The critical noise Pc denotes Shannon's channel capacity, (b) 
K = 2 and L = 2; the ferromagnetic solution and its mirror image are the only minima of the free 
energy over a relatively small noise level (the solid line in the horizontal). The critical point, due to 
dynamical considerations, is the spinodal point Ps where sub-optimal ferromagnetic solutions (F', 
m < 1) emerge. The thermodynamic transition point p^, at which the ferromagnetic solution loses 
its dominance, is below the maximum noise level given by the channel capacity, which implies that 
these codes do not saturate Shannon's bound even if optimally decoded, (c) K = 1; the solid line 
in the horizontal axis represents the range of noise levels where the ferromagnetic state (F) is the 
only minimum of the free energy. The sub-optimal ferromagnetic state (F') appears in the region 
represented by the dashed line. The spinodal point Ps, where F' solution first appears, provides the 
highest noise value in which convergence to the ferromagnetic solution is guaranteed. For higher 
noise levels, the system becomes bistable and an additional unstable solution for the saddle point 
equations necessarily appears. A thermodynamical transition occurs at the noise level pi where 
the state F' becomes dominant. 



31 




FIG. 2. Free energies obtained by solving the analytical equations using Monte-Carlo integra- 
tions for K = 1, R = 1/6 and several values of L. Full lines represent the ferromagnetic free energy 
(FERRO, higher on the right) and the suboptimal ferromagnetic free energy (higher on the left) 
for values of L = 1, 7. The dashed line indicates Shannon's bound and the arrows represent the 
spinodal point values ps for L = 2, 7. The thermodynamic transition is very close, but bellow, 
the channel capacity {pi « 0.261 against pc ~ 0.264 at R= 1/6). 



32 




0.1 0.2 0.3 0.4 0.5 
P 

FIG. 3. Critical code rate as a function of the flip rate p, obtained from numerical solutions 
and the TAP approach (A^ = 10^), and averaged over 10 different initial conditions with error bars 
much smaller than the symbols size, (a) Numerical solutions for if = L = 3, C = 6 and varying input 
bias fs (□) and TAP solutions for both unbiased (+) and biased (O) messages; initial conditions 
were chosen close to the analytical ones. The critical rate is multiplied by the source information 
content to obtain the maximal information transmission rate, which clearly does not go beyond 
i?=3/6 in the case of biased messages; for unbiased patterns H2{fs) = 1- (b) For the unbiased case 
of K = L = 2; initial conditions for the TAP (+) and the numerical solutions {O) were chosen to 
be of almost zero magnetisation, (c) For the case of K = 1, L = 2 and unbiased messages. We 
show numerical solutions of the analytical equations (O) and those obtained by the TAP approach 
(+). The dashed line indicates the performance oi K = L = 2 codes for comparison. Codes with 
K = 1, L = 2 outperform if = L = 2 for code rates R< 1/3. 



33 



0.25 

0.2 

</} 

0.15 
0.1 



"I 1 1 r 

Shannon's bound 



4 5 
L 



FIG. 4. The spinodal point noise level Ps iov K = 1, R = 1/6 and several choices of L. 
Numerical solutions are denoted by circles and TAP decoding solutions [N = 1{)^) by black triangles. 
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TABLES 



TABLE L Comparison between the maximal tolerable noise level for codes based on randomly 
and systematically structured matrices in the case of AT = L = 2; decoding is carried out using 
BP/TAP and the transmission channel used is the BSC. The performance of both matrix structures 
is highly similar. 



Rate R = K/C 


0.6666 


0.5 


0.4 


0.3333 


0.2857 


0.2 


0.1 


Systematic Matrix 


0.0527 


0.0934 


0.1222 


0.1416 


0.1598 


0.1927 


0.2476 




±0.0016 


±0.0019 


±0.0012 


±0.0016 


±0.0007 


±0.0016 


±0.0010 


Random Matrix 


0.0528 


0.0930 


0.1206 


0.1439 


0.1599 


0.1931 


0.2477 




±0.0009 


±0.0019 


±0.0010 


±0.0017 


±0.0010 


±0.0014 


±0.0014 



35 



